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ABSTRACT 


The  research  program  in  seismic  exploration  in  progress  in  the 
Mathematics  Department  of  the  University  of  Denver  is  described.  This 
work  is  identified  here  by  the  term  velocity  inversion.  The  mathematical 
formulations  employed  by  this  group  are  outlined  and  results  of  computer 
implementation  are  depicted.  Ongoing  research  is  also  presented. 


1.  Introduction 

The  purpose  of  this  paper  is  to  describe  the  research  efforts  of  a 
group  of  people  in  the  Mathematics  Department  at  the  University  of  Denver 
and  related  work,  elsewhere.  The  researchers  in  the  group  include  the 
authors.  Professor  F.  G.  Hagin  of  the  Mathematics  Department,  former  gra¬ 
duate  students,  Drs.  R.  D.  Mager,  J.  A.  Armstrong  and  S.  B.  Gray,  present 
graduate  student,  M.  Lahlou  and  programmer,  W.  S.  Grady.  The  objective 
of  this  research  is  to  study  the  problem  of  the  mapping  of  the  interior 
of  the  earth  as  an  inverse  problem  and  to  develop  methods  which  yield 
increasingly  more  accurate  solutions  of  that  inverse  problem.  This 
presentation  was  originally  prepared  as  a  lecture  which  was  presented  at 
the  Seismic  Inversion  Workshop  at  the  50th  Annual  International  Meeting 
of  the  SEG,  Houston,  November  20,  1980. 

The  methods  we  use  are  classical,  employing  perturbation  techniques, 
transform  methods,  asymptotic  and  numerical  analysis  to  arrive  at  com¬ 
puter  algorithms  which  produce  a  mapping  of  the  interior  of  the  earth. 

We  shall  describe,  here,  what  we  mean  by  inversion  in  contrast  to 
migration.  Our  mathematical  development  and  subsequent  computer  imple¬ 
mentation  will  be  presented,  followed  by  a  brief  discussion  of  present 
and  ongoing  research. 


2 


The  mapping  of  the  interior  of  the  earth  from  observations  on  the 
surface  of  the  earth  is  an  inverse  problem.  For  this  type  of  inverse 
pioblem,  the  propagation  of  signals  -  acoustic,  elastic  or  electromag¬ 
netic  -  into  the  earth  is  modelled  by  the  appropriate  equation  or  system 
of  equations  in  which  one  or  more  functions  characterizing  the  interior 
of  the  earth  (soundspeed,  elastic  or  electromagnetic  coefficients)  are 
left  free.  One  or  more  signals  consistent  with  the  model  are  introduced 
at  or  near  the  surface  of  the  earth  in  a  region  of  interest.  The  'irre¬ 
gularities'  of  the  interior  of  the  earth  produce  a  'response'  to  those 
signals.  Observations  of  those  responses  are  recorded. 

The  objective  of  the  inverse  problem  is  to  determine  the  free  coef¬ 
ficients  in  the  modelling  equations  from  knowledge  of  the  input  signal(s) 
and  the  response (s)  and  thereby  'map'  the  interior  of  the  earth. 

This  type  of  inverse  problem  is  known  by  the  fuller  title,  inverse 
scattering  problem.  This  contrasts  with  the  more  familiar  direct 
scattering  problem  in  which  the  parameters  of  the  equation(s)  are  known 
and  the  objective  is  to  determine  the  response  to  the  given  signal 

The  mapping  of  the  interior  of  the  earth  from  observations  of  the 
response  to  a  single  acoustic  source  is  an  inverse  problem  for  the  acous¬ 
tic  wave  equation.  Here,  implicit  in  the  model  -  acoustic  wave  equation 
-  is  a  definition  of  the  word  'mapping'.  At  most,  one  could  hope  to 


characterize  the  interior  of  the  earth  frost  this  model  in  terms  of  its 
density  and  compressibility  or  its  density  and  acoustic  propagation 
speed. 

Under  the  assumption  of  constant  density,  an  approximate  solution  to 
this  inverse  problem  for  the  velocity  was  demonstrated  by  Claerbout 
[1971]  .  Approximate  solutions  for  both  velocity  and  density  in  a  hor¬ 
izontally  stratified  earth  have  been  presented  by  Raz  [1981c]  and  Clayton 
and  Stolt  [1980]  . 

The  mapping  of  the  interior  of  the  earth  from  'local'  observations 
or.  the  surface  of  the  responses  to  an  'ensemble'  of  acoustic  sources  is 
another  type  of  inverse  problem.  Each  experiment  in  the  ensemble  is  per¬ 
formed  separately  and  the  'locality'  of  the  observations  is  limited  by 
the  extent  of  the  receiver  array. 

This  is  the  inverse  problem  which  has  dominated  exploration  geophy¬ 
sics  in  the  recent  past.  Below,  this  problem  will  be  referred  to  by  the 
title,  the  inverse  problem  of  seismic  exploration,  although  we  ack¬ 
nowledge  the  introduction  into  the  repertoire  of  exploration  geophysics 
of  other  inverse  problems  -  static  and  dynamic,  acoustic,  elastic  and 
electromagnetic. 

The  seismic  exploration  problem  has  recently  been  treated  most  suc¬ 
cessfully  by  wave  equation  migration  [Claerbout  and  Doherty,  1972].  This 
technique  treats  the  inverse  problem  of  seismic  exploration  as  a  direct 
problem  in  reversed  time  for  the  acoustic  wave  equation  with  'halved' 


propagation  speed. 


The  following  features  of  wave  equation  migration  are  to  be  noted: 


It.  is 


i.  a  high  frequency  technique, 

ii.  emphasizing  phase,  but 

iii.  neglecting  amplitude. 

Migration 

i.  locates  reflectors,  but 

ii.  does  not  estimate  'reflection  strength*. 

Our  approach  to  the  inverse  problem  of  seismic  exploration  is  iden¬ 
tified  by  the  name,  velocity  inversion  [See,  Cohen  and  Bleistein,  1977, 
1979a].  In  this  approach,  the  inverse  problem  is  modelled  by  one  or  more 
non-linear  integral  equations  deduced  from  the  direct  scattering  model. 
The  integral  equation(s)  are  linearized  by  perturbation  methods  (or,  as 
the  physicists  would  prefer,  by  a  Born-like  approximation).  The  unknowns 
in  these  equations  are  the  free  parameters  of  the  model  equations.  Thus, 
inversion  of  the  integral  equation  leads  directly  to  a  solution  of  the 
inverse  problem. 
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Velocity  inversion,  as  implemented  on  seismic  data,  is 


i.  a  high  frequency  technique, 

ii.  preserving  phase  and 

iii.  preserving  amplitude. 

Therefore,  velocity  inversion 

i.  locates  reflectors  and 

ii.  estimates  reflection  strength. 

The  accuracy  of  the  estimate  of  reflection  strength  is  limited  by 
the  accuracy  of  the  perturbation  approximation.  Thus,  much  current 
research  by  us  and  by  others,  notably  Coen  [1980],  is  focused  on  updating 
the  perturbation,  or  correcting  it,  or  perturbing  abont  more  accurate 
reference  signals. 

It  should  be  further  noted  that  high  frequency  implementation  is  not 
a  constraint  of  the  basic  method,  but  rather  an  acknowledgement  of  the 
realities  of  seismic  data.  For  the  length  scales,  velocities  and  curva¬ 
tures  of  interest,  4  Hz  is,  indeed,  a  'high  frequency'.  For  example,  for 
a  length  scale,  L*=1000  ft.,  a  velocity  of  v=5000  ft/sec,  and  a  frequency, 
f=4  Hz,  the  relevant  parameter. 
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2nfL/v=5, 


is  ' large ' . 

Finally,  we  remark  that  estimation  of  velocity  from  information 
about  reflection  strength  requires  a  further  assumption  about  the  rela¬ 
tion  between  density  and  compressibility.  We  assume  a  constant  density 
medium. 
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For  analysis  of  the  inverse  problem  of  seismic  exploration  by  velo¬ 
city  inversion,  the  propagation  of  acoustic  waves  in  the  interior  of  the 
earth  is  modelled  by  the  homogeneous  acoustic  wave  equation, 

V*D  -  =0,  z  >  0,  (3.1) 

subject  to  the  introduction  of  appropriate  sources  on  the  upper  surface 
[See,  Cohen  and  Bleistein,  1977,  1979a].  It  is  assumed  that  everywhere 
on  the  surface,  a  'backs ca  ttered  response*  is  observed.  That  is,  the 
upward  travelling  signal  is  observed  at  the  source  point.  This  is  the 
mathematical  idealization  of  an  array  of  CDP  stacked  responses. 

The  velocity  v  is  rewritten  as 


-i=^(l+a).  (3.2) 

v*  c1 

An  integral  equation  is  then  derived  for  a  by  a  standard  perturbation 
technique.  That  integral  equation  is 


dx  dy  a(x,y, z)c-*(x,y,z) 


Dj(t, x,y, z, ^ ,n,0)Uj(r-t,5,q,0;x,y, z) 
Us(t,?,n,0;5»il>0)(T-t).  (3.3) 


Here  Uj  is  the  impulse  response  due  to  a  source  at  the  point  (£,q,0) 


at  time,  t=0,  in  the  absence  of  the  perturbation,  a.  U_  is  the  observed 
'backscattered'  wave  at  due  to  the  presence  of  the  perturbation,  a. 

When  the  reference  velocity,  c,  is  assumed  to  be  constant  and  obser¬ 
vations  of  the  backscattered  field  are  made  'everywhere*  on  the  upper 
surface,  then  this  integral  equation  has  the  following  analytical  solu¬ 
tion: 


a(r,y,z)  =  ^J^dnf^dk^dk,  J^dr  jjdt[k,  (r»-ct) 

Us(t»?»’l»0;§,n»0)exp{2i[k1(x-§)+kl(y-q)-k,z]+i«ot}  ];  (3  .4) 

.  r  2  1 

«  =  ctsign  k,][kl+kJ+k, J 

Thus,  the  solution  consists  of  a  multi-fold  Fourier  transform  over 
the  observations  in  and  t,  followed  by  a  transformation  of  transform 

variables  from  kj.kj.oj  ,  to  k1,k2,k,,  via  the  indicated  dispersion  rela¬ 
tion.  This  is  a  full  bandwidth  solution  for  the  perturbation  in  velo¬ 
city,  a. 


This  formula  is  also  quite  similar  to  the  result  of  F-K  migration 
[Stolt,  1978].  Indeed,  the  methods  differ  only  in  the  amplitude  weight¬ 
ing  factor,  here  deduced  by  solving  exactly  an  approximate  integral  equa¬ 
tion  for  the  true  perturbation  in  velocity. 

Thus,  this  solution  could  be  implemented  by  use  of  fast  Fourier 
transform  in  three  variables,  followed  by  an  interpolation  to  obtain  the 
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data  ot>  the  appropriate  grid,  then  followed  by  the  inverse  Fourier 
transform . 

The  direct  calculation  of  a  in  this  way  does  not  exploit  the  partic¬ 
ular  properties  of  the  data  or  of  the  nature  of  the  solution  we  seek. 
Firstly,  as  already  noted,  the  data  is  bandlimited  to  high  frequency. 
Thus,  the  Fourier  transform  of  a,  implicit  in  (3.4)  before  calculation  of 
the  inverse  transform,  is  bandl imited . 

Tn  a  series  of  papers  [Bojarski,  1966;  Mager  and  Bleistein,  1978; 
Armstrong  and  Bleistein,  1978;  Cohen  and  Bleistein,  1979b],  a  theory  was 
developed  to  extract  information  from  a  high  frequency  bandlimited 
Fourier  transform  of  a  piecewise  constant  function.  More  precisely,  it 
is  shown  in  this  series  of  references  ho  ’  to  locate  the  discontinuities 
of  such  a  function  and  how  to  estimate  the  magnitude  of  the  discontinuity 
from  the  Fourier  inversion  of  the  bandlimited  data.  To  obtain  this 
information,  the  Fourier  data  is  processed  to  yield  the  normal  der’vative 
to  the  surface(s)  of  discontinuity.  This  derivative  is  a  bandlimited 
Dirac  delta  function,  which  is  readily  recognizable.  The  ’strength'  of 
the  delta  function  at  its  peak  is  in  known  proportion  to  the  magnitude  of 
the  discontinuity  and  to  the  bandwidth. 

For  the  seismic  inverse  problem,  the  output  of  this  process  is  an 
array  of  delta  functions  which  define  the  'events'  or  boundaries  between 
the  layers  of  the  subsurface.  From  the  peak  values  of  the  output,  the 
reflection  strength  or  velocity  increments  can  be  calculated. 


For  the  integral  in  (3.1),  the  fact  that  only  high  frequency  data  is 


usually  collected  can  be  exploited  to  reduce  the  number  of  integrations 


to  be  performed.  With  fewer  integrations,  it  is  practical  to  develop  an 
algorithm  which  computes  the  output  pointwise.  Such  a  technique  has  the 
advantage  that  the  velocity  need  not  be  kept  constant  for  each  point  of 
processing,  but  can  be  replaced  by  a  'local'  estimate  of  the  rms  velo¬ 
city.  Implementation  of  such  a  variable  reference  velocity  tends  to 
place  the  depth  of  events  of  the  output  more  realistically  and  also 
allows  diffraction  tails  ’at  depth'  to  be  gathered  up  more  completely. 
An  example  of  the  latter  will  be  shown  below. 

It  should  be  noted  that  (3.4)  provides  a  three  dimensional  output 
under  the  assumption  that  a  full  two  dimensional  ensemble  of  observations 
is  made  at  the  upper  surface.  More  typically  in  seismic  exploration,  a 
line  of  data  is  taken.  It  is  then  assumed  that  the  earth  has  only  two 
dimensional  variation  -  along  that  line  and  vertically.  This  is 
equivalent  to  assuming  y- independence  of  a  and  q- independence  of  the 
observations.  This  results  in  the  elimination  of  the  q  and  k2  integra¬ 
tions  in  (3.4)  and  also  the  elimination  of  one  factor  of  n.  It  should  be 
noted  that  this  is  the  result  of  assuming  three  dimensional  propagation 
from  a  point  source  over  an  earth  with  two  dimensional  velocity  varia¬ 
tion.  Thus,  the  velocity  increments  are  still  estimated  by  the  output, 
but  only  to  the  extent  that  the  model  is  consistent  with  reality.  This 
is  in  contrast  to  a  completely  two  dimensional  model  in  which  point 
sources  are  equivalent  to  1  ine  sources  over  the  subsurface.  In  that 
case,  even  if  the  earth  did  have  only  two  dimensional  variation,  the  out¬ 
put  would  not  properly  estimate  velocity  increments,  because  the  sources 
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Hat, tmm 


and  the  nature  of  propagation  in  the  three  dimensional  earth  vas  not 
modelled  correctly. 


I 
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4 . Implementation  of  Velocity  Inversion  on  Synthetic  Data 

A  number  of  synthetic  data  examples  for  the  velocity  inversion  for¬ 
malism  were  presented  in  [Cohen  and  Bleistein,  1979a] .  Those  examples 
confirmed  the  validity  of  the  computer  implementation  of  our  method  on 
data  sets  generated  from  simple  reflectors  and  on  one  example  with  multi¬ 
ple  reflectors,  including  an  interior  lens  shaped  region.  Furthermore, 
it  was  seen  from  those  examples  that  the  perturbation  method  produced 
adequate  estimates  of  velocity  increments  which  were  20%  of  the  reference 
velocity. 

Here,  four  synthetic  examples  will  be  presented.  The  first  two  of 
these  are  examples  with  two  dimensional  variation,  the  last  two  are  exam¬ 
ples  with  three  dimensional  variation. 

Figure  1  is  the  synthetic  timelog  for  the  impulse  response  from  a 
two  dimensional  buried  focus.  The  data  was  generated  from  a  Kirchhoff 
representation  of  the  upward  scattered  wave.  The  doublet-like  behavior 
of  the  response  on  the  lower  curve  of  the  'bowtie'  can  be  seen.  Indeed, 
the  objective  of  this  example  was  to  test  that  velocity  inversion  unrav¬ 
els  this  doublet. 

The  output  of  our  algorithm  is  shown  in  Figure  2.  The  location  of 
the  surface,  as  defined  at  each  point  by  the  peak  of  the  bandlimited 
delta  function,  is  virtually  exact.  For  this  example,  the  velocity  in 
the  upper  medium  was  taken  to  be  5000  ft/sec  and  the  increment  at  the 
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interface  was  taken  to  be  250  ft/sec.  Typical  estimates  of  the  velocity 
increment  taken  from  the  output  were  248  ft/sec  to  251  ft/sec. 

The  timelog  for  the  second  example  is  shown  in  Figure  3.  Here,  the 
reflector  has  a  sharp  discontinuity  and  the  objective  was  to  test  how 
well  the  method  'gathered  up'  the  diffraction  tail. 

The  output  of  our  method  is  shown  in  Figure  4.  Again,  the  location 

of  the  reflector  is  exact.  The  velocity  and  the  velocity  increment  of 

the  synthetic  were  as  in  the  previous  example.  Away  from  the  edge  of  the 
reflector,  the  estimates  of  velocity  increment  were  as  in  the  previous 
example.  At  the  edge  of  the  reflector,  the  estimate  decreased  to  half 
its  true  value,  as  predicted  by  theory.  One  trace  further  to  the  left, 
it  reduced  to  one  tenth  of  the  value  on  the  reflector,  then  to  one  one- 
hundredth  and  then  was  lost  in  the  noise.  Thus,  the  edge  was  adequately 
reproduced,  along  with  an  estimate  of  velocity  increment. 

The  final  two  examples  serve  as  a  test  of  the  computer  program 

developed  to  perform  three  dimensional  velocity  inversion.  The  input  was 
analytically  generated  data  for  the  reflection  from  a  planar  reflector 
and  a  spherical  reflector.  The  output  is  shown  in  Figure  5  (for  the 

planar  reflector)  and  in  Figure  6  (for  the  spherical  reflector). 

For  reasons  of  economy,  the  observation  grid  was  coarsened  and  an 
attendant  coarsening  of  the  output  was  observed.  Nonetheless,  the 
results  do  provide  a  test  of  the  basic  theory  as  implemented  by  our  com¬ 


puter  program  for  three  dimensional  velocity  inversion. 


5.  Implementation  of  Velocity  Inversion  on  Real  Data 

Here,  we  show  the  results  of  applying  velocity  inversion  to  real 
data  sets.  Figure  7  shows  the  timelog  of  a  data  set  provided  to  us  by 
Paul  Stoffa  of  Lamont  Observatory.  A  discussion  of  this  data  set  as 
analyzed  by  wave  equation  migration  can  be  found  in  [Herron,  etal., 
1978].  In  Figure  8,  the  results  of  applying  velocity  inversion  to  this 
data  set  are  depicted.  The  figure  is  normalized  by  the  peak  amplitude  of 
the  data  set  with  no  scaling  of  the  output.  Consequently,  the  picture  is 
completely  dominated  by  the  strong  reflection  at  the  ocean  bottom.  In 
Figure  9,  only  the  part  of  the  depth  section  below  3000  meters  is  plot¬ 
ted.  Here,  a  second  reflector,  sloping  downward  to  the  right,  is  visi¬ 
ble.  This  is  the  reflector  R4  in  the  above  cited  reference. 

The  data  provided,  here,  was  not  true  amplitude  data.  Furthermore, 
it  is  relative  amplitude  data  only  to  the  extent  that  relative  amplitude 
is  preserved  by  stacking  and  other  preprocessing.  Thus,  the  amplitude  of 
the  output  is  relative  amplitude  data,  subject  to  the  same  caveat. 
Nonetheless,  treating  it  as  accurate  relative  amplitude  data,  the  output 
can  be  normalized  with  respect  to  the  assumed  known  response  at  the  ocean 
bottom.  It  was  assumed  that  the  velocity  incremented  at  the  ocean  bottom 
from  1500  meters  to  3000  meters.  On  this  basis,  estimates  were  made  of 
the  velocity  change  at  the  reflector  below  -1625  meters.  The  estimates 
varied  from  -150  m/sec  to  -250  m/sec. 

As  a  basis  of  partial  comparison  of  the  results  of  velocity 
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inversion  and  migration,  an  electrostatic  plot  of  the  type  ased  by 
Lamont,  was  generated  for  ns  by  Stoffa.  That  result  is  shown  along  with 
Stoffa's  migrated  section,  in  Figure  10.  The  similarity  of  the  output  is 
quite  apparent.  That  is,  inversion  provides  as  qualitatively  ’clean'  a 
depiction  of  the  depth  section,  while  also  providing  estimates  of  velo¬ 
city  changes. 


Figures  11  and  12  are  the  left  and  right  halves,  respectively,  of  a 
timelog  provided  to  us  by  Marathon  Oil  Company.  Figures  13  and  14  depict 
the  results  of  velocity  inversion.  In  Figure  13,  a  long  diffraction  tail 
from  the  edge  of  the  reflector  in  Figure  14  can  be  seen.  The  reason  that 
this  diffraction  tail  ’survived*  inversion  became  apparent  in  our  discus¬ 
sions  with  the  geophysical  research  group  at  Marathon.  The  constant 
reference  velocity  used  in  this  first  'pass'  at  the  data  was  5000  ft/sec. 
While  this  was  a  good  estimate  near  the  surface,  8000  ft/sec  was  a  better 
estimate  at  the  depth  of  the  discontinuous  reflector.  Figures  15  and  16 
depict  a  reprocessing  of  the  data  below  4900  ft.  at  this  corrected  speed. 
Now,  it  can  be  seen  that  the  major  diffraction  tail  has  been  'gathered 
up'  .  There  is,  unfortunately,  now  an  ’overmigration'  phenomenon  apparent 
with  respect  to  other  reflectors  at  lower  depth.  A  full  analysis  of  this 
data  would  require  a  number  of  different  reference  velocities  in  dif¬ 
ferent  regions.  We  did  not  have  the  resources  to  carry  this  out.  How¬ 
ever,  the  output  presented  does  demonstrate  the  practicality  of  using 
different  reference  velocities  at  different  depths. 


It  should  be  noted  that  in  this  'atter  processing  step  at  the  new 


velocity,  it  was  not  necessary  to  reprocess  the  portion  of  the  section 


■  -at*,-'* ; . 


above  4900  ft,  but  only  to  process  the  data  in  the  region  of  interest. 


This  would  be  less  practical  in  a  multi-fold  Fourier  inversion,  since 
each  separate  reference  velocity  would  require  a  different  interpolation 
grid  when  transforming  from  k1(k,,w,  to  kj.kj.k,. 

While  this  extension  will  help  clean  up  diffraction  tails,  it  will 
not  correct  for  refractions,  since  the  closed  form  solution  asstnes  an 
incident  wave  which  propagates  on  straight  rays. 

The  number  of  data  points  and  processing  times  for  this  output  are 
listed  in  the  table  below: 


Process 

ing  Time  in 

seconds  on  a 

CDC  Cyber  76 

#  of  input 

#  of  output 

Preproce  ss 

Process 

Total 

Source 

data  points 

data  points 

time(sec) 

time(secs) 

time 

Lamont 

480,000 

60,000 

31 

201 

232 

Marathon 

700,000 

230,000 

117 

583 

700 
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6..  Recent  Development s  and  Future  Research 

The  objectives  of  recent  research  efforts  have  been  to  develop 
methods  to  invert  the  integral  equation  (6.2)  when  the  reference  speed, 
c,  is  not  constant  and  to  more  adequately  account  for  the  non-linearity 
of  the  underlying  inverse  problem. 

6.1  Higher  Order  Accuracy  in  the  One-Dimensional  Problem 

An  important  line  of  research  leading  to  higher  order  accuracy  in 
the  one-dimensional  velocity  inversion  problem  was  initiated  indepen¬ 
dently  and  nearly  simultaneously  by  Gray  [1980]  and  Raz  [1981c].  While 
their  approach  differs,  the  basic  result  obtained  is  nearly  identical. 

We  shall  describe  Gray's  approach.  He  introduces  the  travel  time  as 
an  independent  variable.  The  resulting  equation  admits  an  incident  wave 
with  constant  reference  speed  and  a  perturbation  proportional  to  the  log¬ 
arithmic  derivative  of  the  true  propagation  speed.  The  integral  equation 
for  this  perturbation  reduces  to  an  equation  relating  its  Fourier 
transform  to  the  observed  backsca ttered  data.  Fourier  inversion  and 
exponentiation  then  produces  the  velocity  a_s  &  function  of  travel  time . 
Since  the  true  depth  is  an  integral  of  this  velocity  with  respect  to 
travel  time,  the  velocity  is  given  implicitly  in  terms  of  depth. 


Gray  shows  that  this  result  has  higher  order  accuracy  in  the  pertur¬ 
bation,  a.  Gray  and  Hagin  [1980]  have  shown  that  an  iteration  scheme 
based  on  this  formulation  converges  even  when  the  velocity  has  jump 


discontinuities  (albeit,  not  too  large),  as  in  a  layered  earth.  In  con¬ 
trast,  Prosser  [1980]  has  shown  that  an  iteration  scheme  based  on  the 
previous  perturbation  technique  will  converge  for  sufficiently  smooth 
velocities.  However,  in  private  communications,  he  has  indicated  pessi¬ 
mism  about  demonstrating  convergence  of  that  scheme  for  the  layered 
medium  case,  which  arises  in  the  seismic  exploration  problem. 

Figure  17  is  taken  from  Gray,  Bleistein  and  Cohen  [1980] .  It  demon¬ 
strates  the  increased  accuracy  of  this  method,  both  in  locating  the 
discontinuities  and  estimating  their  magnitude.  Figure  18,  taken  from 
Ilagin  [1980]  demonstrates  the  increased  accuracy  of  second  iterates  in 
this  method. 


Raz  [1979,  1981a]  has  provided  some  extensions  of  this  work  to  the 
three-dimensional  stratified  case.  He  has  shown  how  observations  at 
offset  can  be  used  to  invert  a  tilted  stratified  earth  or  could  be  used 
to  separate  density  and  compressibility.  An  alternative,  but  quite  simi¬ 
lar  approach,  to  the  latter  problem  was  also  presented  at  the  SEG  meeting 
in  Houston  by  Stolt  and  Jacobs  [1980]. 


Our  research  group  is  presently  investigating  the  extension  of  this 
method  to  the  three  dimensional  seismic  inverse  problem. 

6.2  A  Wave  Equation  for  the  Propagation  of  the  Ensemble  of  Backscattered 
Signals 


I 

I 


Research  on  an  alternative  approach  to  the  three  dimensional  velo¬ 
city  inversion  problem  has  recently  been  initiated.  The  objective  of 
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this  approach  is  to  derive  an  equation  for  the  propagation  of  the  ensem¬ 
ble  of  backscattered  signals.  That  equation  should  be  sufficiently  accu¬ 
rate  to  preserve  both  phase  and  amplitude  of  the  ensemble.  Clearly,  this 
is  the  objective  of  wave  equation  migration.  However,  that  method  is 
based  on  travel  time  arguments  alone  and  thus  accurately  produces  only 
phase  information,  except  at  the  first  reflector.  This  will  be  discussed 
in  further  detail,  below. 

Our  method  is  based  on  analysis  of  the  Kirchhoff  integral  represen¬ 
tation  of  the  backscattered  impulse  response.  For  a  single  layer,  we 
consider  the  geometrical  configuration  depicted  in  Figure  19.  By  using 
geometrical  optics  approximations  in  the  Kirchhoff  integral  representa¬ 
tion  of  the  solution,  one  can  derive  the  following  integral  representa¬ 
tion  for  the  backscattered  wave  (See,  for  example,  Cohen  and  Bleistein 
1979b): 


us(x,w) 


i<o  f 
8n*c  Js 


exp{2iu>/c) 


dS. 


(6.1) 


Here,  2  is  the  unit  upward  normal  vector.  T  is  a  unit  vector  from  the 

surface,  S,  to  the  observation  point,  u  denotes  the  Fourier  time 

o 

transform  of  Uc  and  R  is  the  geometrical  optics  reflection  coefficient, 
2.7  -  (c*/c11  -  1  +  (2.7)*)l/* 

R(x,£)  =  - jy  .  (6.2) 

n«7  +  (c*/Cl*  -  1  +  (2.?)*)  '* 


By  applying  the  wave  operator  to  (6.1)  and  retaining  only  terms  to 
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two 


orders  in  <o  (i.e.,  terms  multiplied  by  and  u  ),  the  following 


partial  differential  equation  can  be  derived  for  U^: 


[5.  *,»(.,). 


Here,  6  denotes  the  Dirac  delta  function*  a  denotes  normal  distance 

n 

from  the  reflector  to  the  observation  point.  R  denotes  the  'normal' 

n 

reflection  coefficient: 


ct  ~  c 

R  =  - 

n  Cj  +  c 


(6.4) 


Boundary  data  for  u^  is  also  prescribed  at  z=0. 

Equation  (6.3)  has  the  interpretation  that  the  ensemble  of  back- 
scatters  is  the  response  to  a  source  distributed  over  the  reflecting  sur¬ 
face.  Furthermore,  that  source  is  proportional  to  the  reflecting 
strength.  Thus,  equation  (6.3)  provides  a  'quantification'  of  the 
'explosive  reflector'  model  of  backsca ttering  and  gives  the  intensity  of 
the  source  to  leading  order  in  u). 

The  inverse  problem  for  the  location  of  the  reflector  and  the 
reflection  strength  may  now  be  stated  as  follows: 


Since  the  'support'  of  the  delta  function  is  on  the  scattering 
surface,  this  distance  need  only  be  single  valued  and  well  defined 
near  that  surface,  which  it  will  be  for  a  sufficiently  'well 
behaved'  surface. 


let  ug  be  a  solution  of  (6.3)  with  unknown  source.  Given  the  boun¬ 
dary  data  for  ug,  find  the  source. 

Unfortunately,  inverse  source  problems  are  known  to  have  non-unique 
solutions.  7t  is  our  belief,  however,  that  in  the  class  of  distribu¬ 
tional  solutions  with  uniform  strength  over  a  surface,  the  solution  is 
unique.  However,  no  proof  has  been  developed,  as  yet. 

This  result  can  also  be  stated  in  time  domain  with  an  alternative 
interpretation.  To  deduce  this  result,  we  introduce  the  function 

0US 

v(x,<o)  =  -i  (6.5) 

and  its  inverse  Fourier  transform, 

V(x,t)  =  tUs(x,t).  (6.6) 
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This  function  is  a  solution  of  the  following  initial-boundary  value 


problem  deduced  from  the  frequency  domain  problem: 


w  -  erg? 


0, 


V(x.O) 


R  6 (o 
n 

4jt 


Hi 


£V(x,0)  =  0 

at 


(6.7) 


(6.8) 


subject  to  an  appropriate  boundary  condition  on  the  surface  z=0. 

The  inverse  problem  for  the  reflector  and  the  reflection  strength 
may  now  be  stated  as  follows: 


Let  V  be  a  solution  of  the  wave  equation  with  propagation  speed  c/2. 
Given  the  observations  at  the  upper  surface,  z=0,  for  all  time  and 
given  the  second  initial  condition  in  (6.8),  determine  the  initial 
value  of  V. 


This  is  exactly  the  problem  which  is  solved  by  wave  equation  migra¬ 
tion.  We  conclude,  therefore,  that  at.  .the  'first  reflector',  wave  equa¬ 
tion  migration  not  only  locates  the  reflector  but  also,  for  true  ampli¬ 
tude  data,  accurately  estimates  the  reflection  strength. 

It  is  possible  to  obtain  an  integral  representation  of  a  solution  to 
this  problem,  namely. 
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(6.9) 


V<£,0)  =  [  dxdy  j  dt  V(x,y,0,t)  Y>0>L>  *)  . 

Jz=0  Jo  5z 

Here,  G  is  the  half  space  Green's  function  for  the  wave  equation  with 
source  point,  (x,y,0)  and  observation  point,  (£)  and  zero  boundary  value 
at  z=0 . 


In  order  to  generalize  this  result  to  a  layered  earth,  it  is  neces¬ 
sary  to  analyze  the  backsca ttered  field  from  a  layer  below  the  first.  We 
start  from  the  same  Kirchhoff  integral  representation  of  u  However, 
the  function, 

[4nr)~1exp{i(or/c} , 

whose  square  appears  in  the  integral  in  (6.1),  must  now  be  replaced  by 
the  Green's  function  for  an  arbitrary  medium.  This  can  be  simplified 
somewhat  if  we  are  willing  to  accept  only  the  primary  downward  propagat¬ 
ing  part  of  the  Green's  function  and  then  further  content  ourselves  with 
a  characterization  of  that  function  which  has  accurate  phase  and  leading 
order  amplitude,  only.  Thus,  each  reflector  will  be  treated  as  if  it 
were  the  only  reflector;  i.e.,  multiple  reflections  are  neglected  and,  as 
above,  the  representation  reproduces  the  phase  and  amplitude  to  leading 
order  in  w. 

This  asymptotic  Green's  function  has  the  form. 
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G  ~  A  exp{i<o^} 


(6.10) 


where  f  and  A  are  solutions  of  the  eikonal  equation  and  transport  equa¬ 
tion.  respectively, 

(V*)1  =  c-*i  2V*.VA  +  AV»?  =  0.  (6.11) 

Furthermore,  in  order  that  this  be  asymptotically  the  Green's  function, 
it  is  required  that  A  behave  as  [4nr]-1  when  r,  the  distance  between 
source  and  observation  point,  approaches  zero  and  that  f  be  zero  in  that 
limit.  In  this  case,  in  analogy  with  (6.1),  the  backscattered  field  has 
the  integral  representation, 

Ug(x,(i))  ~  2iwj^R  3 •  Vp  A1  exp(2iwp)  dS.  (6.12) 

Here,  the  reflection  coefficient  is  as  in  (6.2)  except  that  7  must  be 
replaced  by  cVp. 

Applying  the  wave  operator  to  (6.12)  yields  the  following  equation: 

[v  *  -  -  ^1  *.»<•.> 

+  R(n-V^)  A  V^-(AV^)  exp(2iw#}  dS.  (6.13) 

Were  the  second  line  absent  from  this  equation,  then  the  downward  propa- 
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gation  of  the  ensemble  of  backscatters  from  any  layer  would  be  exactly 
the  same  as  those  from  the  first  layer.  In  this  case,  wave  equation 
migration  would  produce  true  amplitude  at  each  reflector  and  would  then. 


indeed,  be  inversion.  The  second  line  in  this  equation  therefore 
represents  the  extent  to  which  downward  propagation  according  to  the  wave 
equation  with  speed  (c/2)  of  the  ensemble  of  backsca ttered  signals  devi¬ 
ates  from  producing  a  'true'  ensemble  of  backscattered  signals  'at 
depth'.  It  should  be  noted  that  this  second  line  has  a  multiplier  of  u 
(rather  than  <■>*)  and  thus  asymptotically  affects  the  leading  order  ampli¬ 
tude  and  not  the  phase  of  the  downward  continued  ensemble.  Thus, 
neglecting  this  term  will  cause  errors  in  amplitude  alone  and  not  in 
phase.  This  is  consistent  with  the  empirical  results  demonstrating  that 
wave  equation  migration  (i.e.,  solving  the  equation  represented  by  the 
first  line  only,  in  (6.13))  locates  reflectors  accurately. 

Research  on  this  full  propagation  equation  is  now  being  carried  out 
by  our  research  group. 
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2.  Conclnsions 


A  summary  of  a  research  program  on  a  class  of  methods  collectively 
referred  to  as  velocity  inversion  has  been  presented.  The  mathematical 
formulation  and  computer  implementation  have  been  described  id  new  and 
incomplete  research  in  progress  has  been  discussed. 

An  inspection  of  the  program  for  the  47th  Annual  International  Meet¬ 
ing  of  the  SEG  in  Houston  in  1976  will  reveal  a  paucity  of  papers  on  the 
subject  of  inversion.  In  contrast,  the  50th  Annual  Meeting  in  1980  had 
three  sessions  and  one  workshop  on  inversion.  Having  presented  one  of 
the  papers  on  inversion  in  that  47th  Meeting,  we  view  this  as  a  favorable 
and  gratifying  trend. 
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FIGURE  CAPTIONS 


Fig.  1  Synthetic  tiaelog  from  a  two  dimensional  buried  focus.  Hor¬ 
izontal  spacing,  100  ft.;  bandwidth,  12-36  hz.  Doublet 
behavior  on  lower  half  of  bowtie  is  evident. 

Fig.  2  Output  of  the  velocity  inversion  algorithm  for  the  line  of  data 
of  Fig.  1.  Horizontal  spacing  between  traces,  100  ft.  Doublet 
in  timelog  has  been  replaced  by  impulse  on  the  output  section. 

Fig.  3  Synthetic  timelog  for  a  half-plane  reflector.  Horizontal  spac¬ 
ing,  100  ft.;  bandwidth,  12-36  hz.  Diffraction  tail  to  left  of 
reflector  is  evident. 

Fig.  4  The  output  of  the  velocity  inversion  algorithm  for  the  dataset 
of  Fig.  3.  The  diffraction  tail  has  been  'gathered'  by  the 
algorithm. 

Fig.  S  Output  of  the  three  dimensional  velocity  inversion  algorithm 

for  a  planar  array  of  synthetic  data  over  a  tilted  planar 
reflector.  Horizontal  spacing  is  100  ft.  in  each  transverse 
direction. 

Fig.  6  Output  of  the  three  dimensional  velocity  inversion  algorithm 

for  a  planar  array  of  synthetic  data  over  a  spherical  reflec¬ 
tor. 

Fig.  7  Timelog  of  a  dataset  gathered  by  Lamont-Doher ty  Observatory 

over  the  east  pacific  rise.  Horizontal  spacing,  50  m; 
bandwidth,  8-31hz.  Electrostatic  plot  provided  by  Paul  Stoffa. 

Fig.  8  Output  of  the  two  dimensional  velocity  inversion  algorithm 

applied  to  the  dataset  of  Fig.  7. 

Fig.  9  Replotting  of  the  section  of  Fig.  8  below  3000m. 

Fig.  10  Electrostatic  plot  of  output  of  Fig.  8  and  electrostatic  plot 
of  Stoffa’s  wave  equation  migration  algorithm  applied  to  the 
same  dataset. 

Fig.  11  Left  half  of  timelog  of  a  dataset  provided  by  Marathon  Oil 

Company.  Horizontal  spacing,  82  ft.;  bandwidth,  5-37hz. 

Fig.  12  Right  half  of  dataset  connected  to  left  half  in  Fig.  10. 

Fig.  13  Left  half  output  of  velocity  inversion  algorithm.  Note  long 

diffraction  tail  below  4900  ft. 
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Fig.  14  Right  half  output  of  velocity  inversion  algorithm  adjacent  to 

left  half  of  Fig.  12. 

Fig.  15  Left  half  output  of  velocity  inversion  below  4900  ft.  with  new 
reference  velocity,  9000ft/sec. 

Fig.  16  Right  half  output  of  velocity  inversion  below  4900  ft.  with 

new  reference  velocity,  9000  ft/sec. 

Fig.  17  Comparison  of  first  and  second  order  one-dimensional  velocity 
inversion  algorithm.  Both  amplitude  and  phase  (location  of 
second  discontinuity)  are  more  accurate  with  second  order 
method. 

Fig.  18  Comparison  of  first  and  second  iterates  of  Gray's  second  order 
method. 

Fig.  19  The  geometry  for  analysis  of  the  propagation  of  an  ensemble  of 
backscatters  from  a  single  layer. 
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